#delimit;

****;
version 9.0;

clear;
clear matrix;
estimates clear;
eststo clear;
set more off;
set mem 400m;
set seed 365476247;

****************************************;
*Did not follow;

cap prog drop runme;
prog def runme;

local hypothesis = `0';
tempfile main bootsave;

use "$da/Crimes-Against-Morality_attrition-data.dta", clear;

di ;

areg attrit_sti i_sept_any1 sept_any1 closing, cluster(wsid)
				 abs(wsid);
global coeff2 = "i_sept_any1";
			  
global mainbeta = _b[$coeff2 ];
global maint = (_b[$coeff2 ]-`hypothesis')/_se[$coeff2 ];
predict epshat, resid;
predict yhat, xb;

*also general the "impose the null hypothesis" yhat and residual;
qui gen temp_y = $var - $coeff2 *`hypothesis';
areg temp_y sept_any1 closing, abs(wsid);
predict epshat_imposed, resid;
predict yhat_imposed, xb;
qui replace yhat_imposed = yhat_imposed + $coeff2 *`hypothesis';

sort wsid el closing;
qui save `main', replace;

*getting a count of then number of clusters in the sample;
qui by wsid: keep if _n==1;
qui summ;
global numclust = r(N);

cap erase `bootsave';
qui postfile bskeep beta_np t_np t_wild
	using `bootsave', replace;

forvalues b=1/$bootreps { ;

qui use `main', replace;
qui by wsid: gen tempbs = uniform();
qui by wsid: gen pos = (tempbs[1] <.5);

qui gen wildresid = epshat_imposed *(2*pos - 1);
qui gen wildy = yhat_imposed + wildresid;
qui areg wildy i_sept_any1 sept_any1 closing, cluster(wsid)
				 abs(wsid);
local bst_wild =(_b[$coeff2 ]-`hypothesis')/_se[$coeff2 ];

*nonparametric bootstrap;
bsample, cluster(wsid) idcluster(newsid);
qui areg attrit_sti i_sept_any1 sept_any1 closing, cluster(newsid) abs(newsid);
local bsbeta = _b[$coeff2 ];
local bst = (_b[$coeff2 ]-$mainbeta )/_se[$coeff2 ];

post bskeep (`bsbeta') (`bst') (`bst_wild');
};
*;
qui postclose bskeep;

qui drop _all;
qui set obs 1;
qui gen t_np = $maint;
qui gen t_wild = $maint;
qui append using `bootsave';

qui gen n = .;

foreach stat in t_np t_wild {;
qui summ `stat';
local bign = r(N);
sort `stat';
qui replace n = _n;
qui summ n if abs(`stat'-$maint )< 0.000001;
local myp = r(mean) / `bign';
global pctile_`stat' = 2*min(`myp', (1-`myp'));
};
*;

global mainp = norm($maint);
global pctile_main = 2*min($mainp, (1-$mainp));

local myfmt = "%7.5f";

di;
di "Number BS Reps = $bootreps, Null hypothesis = `hypothesis'";
display "Main beta" _column(13) "main T" _column(22) "Main %le" 
		_column(33) "NP BS %le" _column(44) "wild %le";
di %6.3f $mainbeta _column(13) %6.3f $maint _column(23) 
		`myfmt' $pctile_main _column(34) `myfmt' $pctile_t_np 
		_column(44) `myfmt' $pctile_t_wild ;
		
end;

global bootreps = 999;
runme 0;

